Phase space flows for non-Hamiltonian systems with constraints 



in 
o 
o 

(N 

< 

00 



Alessandro Sergi * 

Dipartimento di Fisica, Sezione Fisica Teorica, Universitd degli Studi di Messina, 
Contrada Papardo CP. 50-98166 Messina, Italy 

In this paper, non-Hamiltonian systems with holonomic constraints are treated by a general- 
ization of Dirac's formalism. Non-Hamiltonian phase space flows can be described by generalized 
antisymmetric brackets or by general Liouville operators which cannot be derived from brackets. 
Both situations are treated. In the first case, a Nose-Dirac bracket is introduced as an example. In 
the second one, Dirac's recipe for projecting out constrained variables from time translation oper- 
ators is generalized and then applied to non-Hamiltonian linear response. Dirac's formalism avoids 
spurious terms in the response function of constrained systems. However, corrections coming from 
phase space measure must be considered for general perturbations. 
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I. INTRODUCTION 



Constrained systems are ubiquitous in theory and com- 
putation and formalisms for their treatment are still 
developed 0, 0, Q. Some time ago Dirac showed 
how to formulate generalized Hamiltonian phase space 
flows which automatically satisfy certain class of con- 
straints 0, IE IE 0- These constraints, which are 
called second class, are specified by a non-zero Pois- 
son bracket with the Hamiltonian. Dirac's investigation 
aimed at finding a quantization procedure for relativistic 
fields which have constraints arising from their Lorentz 
or gauge symmetries. Recently it has been shown 3j 
how Dirac's scheme can be also applied to non-relativistic 
systems, such as those addressed by classical molecular 
dynamics simulations in condensed matter, where con- 
straints are used to describe the topology of molecules 
or rare events |E]- In particular it has been shown how 
Dirac's formalism can be subsumed by means of a gener- 
alized bracket introduced in Refs. [IE 113 • 

The concern of this paper is to generalize the approach 
of 3| to cases where the dynamics is non-Hamiltonian 10, 
IIE lia IJJ, IM |151 ll^ 13 • This is useful since compu- 
tational schemes adopt constraints and non-Hamiltonian 
dynamics at the same time, with the latter implement- 
ing both specific thermodynamic conditions, i.e. constant 
temperature (18,. 19. 20.. .21.]. and nonequilibrium pertur- 
bations not derivable from a Hamiltonian |23 . 123 . 124 . l25l | . 
Building on the results given in Ref. 1^, the discussion 
will be limited to systems with holonomic constraints. 
Nonholonomic constraints (such as those involved in the 
formulation of the isokinetic ensemble 0,0) will not be 
addressed in this paper. Two kinds of non-Hamiltonian 
dynamics will be considered. The first is based on the 
generalized brackets introduced in 0, ^ while the sec- 
ond case arises from Liouville operators which cannot 
be expressed by means of antisymmetric brackets. For 
the first case, the antisymmetric structure of the gen- 
eralized bracket will be used to combine Dirac's the- 
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ory with non-Hamiltonian Nose-Hoover dynamics (more 
general dynamics, such as Nose- Hoover chains |2E|, do 
not introduce any conceptual difference). As an appli- 
cation of the Nose-Dirac phase space flow, the unbiasing 
factor, arising when holonomic constraints are used to 
study rare events 0, is re-derived. The second type of 
non-Hamiltonian dynamics requires a generalization of 
Dirac's scheme in order to project out the constrained 
degrees of freedom from any arbitraty Liouville operator. 
Linear response theory will be reviewed and some fine 
points, which are relevant for analyzing dynamics with 
holonomic constraints (such as in the case of molecular 
systems 26, 27, 28, 29j), will be discussed. In particular, 
it will be shown that correction terms, stemming from 
phase space measure, appear in the response function for 
general forms of perturbations. 

The paper is organized as follows. In Sec. In] Dirac's 
Hamiltonian formalism is briefly reviewed. In Sec. IIIII a 
unified bracket for Nose thermostated dynamics and con- 
straints, producing a non-Hamiltonian Nose-Dirac phase 
space flow, is introduced. As an illustration of the for- 
malism, the Nose-Dirac flow is applied in Appendix ^ 
to the discussion of rare events sampling. In Sec. IIVI 
Dirac's recipe, for projecting out the spurious dynam- 
ics of constrained variables, is first generalized to ar- 
bitrary time-translation operators and then applied to 
non-Hamiltonian Liouville operators which cannot be de- 
rived from brackets. Linear response theory is briefly re- 
viewed, discussing how Dirac's prescription avoids fake 
terms coming from constraints. Nevertheless, it is shown 
that corrections terms in the response function, origi- 
nating from the constrained phase space measure, may 
appear in the general case. 



II. HAMILTONIAN FORMALISM FOR 
SYSTEMS WITH HOLONOMIC CONSTRAINTS 



Consider a system with a conserved energy 7io(x), 
where x — (r,p) denotes the phase space point. To 
formulate phase space equations of motion in the pres- 
ence of mechanical constraints one can follow Dirac's ap- 
proach 0,11 |E IE Si] 

Together with the n constraints in 
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configuration space cr(r) = 0, one has to consider an ad- 
ditional number n of phase space constraints &{r,r) — 0. 
It is useful to let x — (""i denote the entire set of 
2n phase space constraints. The n constraints & — 
are redundant but necessary to set up a phase space pic- 
ture of the dynamics. Following the convention (due to 
Dirac of evaluating derivatives first and imposing 

constraint relations after, these constraints will disappear 
from the equations of motion and will not contribute to 
the phase space measure. 

The equations of motion with constraints may be writ- 
ten as 



2N 



(1) 



where 2N is phase space dimension and 3^ is an anti- 
symmetric tensor defined by 



2N 2n 



B?M = - E E ^-if '-g-Bf, , (2) 

k,l a,p 



with 



B 



1 
1 



(3) 



usually called the symplectic matrix |3Clj| . In order to 
arrive at Eq. Q one has to define 



(4) 



given in terms of Poisson bracket of the 2n phase space 
constraints, and the inverse matrix (C""'^)^^ (a,/? = 
1, 2n), written explicitly in block form as 



Z 



z-^ 



where the matrices 
and 



N 



nil 



(5) 



(6) 



N 



Ta/S = E 



Pi 



i,k=l 



rriimk 



^ki^^aVk^f) - VkCr»'^kt'^l3 (7) 



(with a, P — 1, . . . , n) have been defined. 



The matrix can be written explicitly in block form 



(8) 



where 



N n 
k,l — l Q:,/3— 1 

{i,j = l,...,N) 



(9) 



and 



N n 



daa 



k,l=l Q,/3=l 

(i,j = l,...,7V) 



(10) 



Substituting B^ into Eq. JQ) and taking into account 
the fact that & — 0, one obtains the equations of mo- 
tion 121 



rrii 



where 



Pi = Fi + ViO- ■ A(r,p) , 



A(r,p) = Z-i.{^,Ho} 



(11) 
(12) 

(13) 



Equations (|ll|) - (|12|l have a phase space compressibil- 

ity sua 



«c = -^ln||Z| 



and distribution function 0,01 

^ sino)six)\\z\\ 



(14) 



(15) 



where S{x) — Yia ^{^a)S{(Ta)- It is worth to remark that 
Eqs. iQJ with the tensor in Eq. and their explicit 
form (|ll|l - p2|) . can be regarded as Hamiltonian since the 
associated generalized bracket 



N 



, ,N ■^-^ da db 



(16) 



where a and b are arbitrary phase space functions, satis- 
fies the Jacobi relation jj, . The generalized bracket 
in Eq. (|16|l has the property of leaving invariant, by con- 
struction, any function of the constraints. 



III. NOSE-DIRAC PHASE SPACE FLOW 

Starting from the structure of either or the as- 
sociated generalized bracket in (|16|) . with B^ given by 
Eq. it is very easy to define non-Hamiltonian equa- 
tions of motion. It suffices to substitute the tensor B'^ in 
Eq. jSJ with a more general antisymmetric tensor B{x) so 
that the Jacobi identity is no longer satisfied. When one 
tries to apply this program to extended system dynam- 
ics, as in the case of Nose-Hoover dynamics, problems 



3 



are encountered since the constraints are usually defined 
only onto a subspace of the extended phase space vari- 
ables. This straightforward approach would make the 
generalized bracket identically zero. One can bypass this 
problem by exploiting the block structure of as given 
in 0. To this end, consider Nose extended phase space 



with coordinates x 
Hamiltonian 



(r, 77, p,Pri), and introduce the Nose 



N 



7Y^ 



i 



2mr 



gkBTrj 



If one defines the antisymmetric matrix 



B 






-1 + A 




1 A 

1 

A p 

1 p 



(17) 
(18) 

(19) 



and uses this matrix in place of B^ either in or in (|16|l . 
then, through the Nose Hamiltonian in Eq. p8(l , one ob- 
tains the desired equations of motion. 





p_ 


(20) 




rrii 




= F,-t-V*T-A(r,p)-p,^ 


(21) 




rrin 


f] 


_ Prj_ 


(22) 




rrir, 




Pv 




(23) 



with F,^ = J2^P^/r 



gksT. 



The phase space flow defined via B in Eq. H19|l con- 
serves the Hamiltonian and any function of the con- 
straints. The equations of motion (Eqs. (|20() - (|23|l ) have 
a compressibility 



^ dB^j on 

^ dxi dxj 



KN 



(24) 



where is given in Eq. (|14|l . The Nose compressibility is 
kn = 37V t) so that the total compressibility of Eqs H2()(l - 
(1^ is 



dln|Z| 
Jt 



dHr 
dt 



(25) 



The primitive function of the compressibility is w{x) = 
— ln|Z| -I- /37ir so that the distribution function in the 
extended phase space is 



/9ND(r,P,?y,P,,) = 5{'H^)5{(T)5{cT)\Z\e 



(26) 



One can easily prove that Eq. 126() provides the distribu- 
tion of a canonical ensemble with constraints. Integrating 
on rj one has 



drj 



9 



(27) 



The constant can be absorbed in the normalization and 
the Gaussian integration on can be easily performed 
so that one obtains 



Pc(r,p) = 5{x)\Z\e 



(28) 



where Tio — X^i Pf /2wi -f ^(r) is the Hamiltonian of the 
physical degrees of freedom. As an example, Nose-Dirac 
flow will be applied in Appendix to the sampling of 
rare events, and it will be shown how to re-derive the 
unbiasing factor first introduced in Ref. |^. 



IV. GENERAL NON-HAMILTONIAN 
DYNAMICS AND CONSTRAINTS 

In the case of general non-Hamiltonian dynamics one 
cannot derive the generator of time translation from gen- 
eralized brackets. Instead one is led to consider a Liou- 
ville operator of the form 

zLp = P}) • 7^ + D({r, p}) • ^ (29) 



defining the time evolution of any arbitrary phase space 
variable a({r, p}) through d = iL^a. The phase space 
incompressibility condition is usually adopted jl^ 



dC 9D 







(30) 



and for simplicity the same will be done here. In molec- 
ular dynamics applications, Liouville operators, having 
the same form as that in Eq. (|29|l . are used to introduce 
time dependent perturbations by means of operators of 
the form 



iLi{t) = T{t)iLp 



(31) 



The unperturbed systems is usually subjected to the ac- 
tion of an operator iLo which is instead derivable from 
some (generalized or Poisson) bracket with the Hamil- 
tonian. Accordingly, the total dynamics is defined via 
the operator iL{t) = iLo + iLi[t). In the presence 
of holonomic constraints, for example describing rigid 
molecules, the formalism of Ref. for the Hamiltonian 
case, or of the previous section, for non-Hamiltonian dy- 
namics, can be used to define an operator iL^ having the 
constraints as conserved quantities. The problem is that 
iL-p and iLj{t) as such do not preserve the constraints 
and could lead to spurious term in the linear response 
derivation, as it will be shown in the following. This fea- 
ture of the formalism is not desiderable since, in actual 
molecular dynamics calculations, the algorithms enforc- 
ing the constraints is used in t he p resence of the per- 
turbation determined by iLi{t) Hi IH IH lli so that 
this perturbation does not violate the constraints in prac- 
tice. The conclusion is that, in order to set up a correct 
formalism, one must project out the dynamics that iLp 
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and iLj{t) spuriously impose on the constraints. To this 
aim, by using a simple extension of Dirac's recipe to gen- 
eral non-Hamiltonian Liouville operators, one can define 
iLfit) as follows 

iLf{t)a = iLi{t)a 

- ^{a.,Xa}{C'~^)ai3{iLi{i)X(3) 

af3 

— J-{t) [iL-pa 

= T{t)iL^a , (32) 

where a is an arbitrary phase space variable, {a,Xa} is 
the Poisson bracket and C is derived by means of Eqs. Q 
and ||SJ). As in the Hamiltonian case Q or in that of the 
non-Hamiltonian bracket, iLf{t) preserves any function 
of the constraints, by construction. Equation H32|) gener- 
alizes Dirac's theory to arbitrary non-Hamiltonian phase 
space flows. 

By means of iLf (t) one can set up the correct for- 
malism for the linear response of systems with holonomic 
constraints subject to a non-Hamiltonian time-dependent 
perturbation. For simplicity, the case in which the un- 
perturbed dynamics of the constrained system is Hamil- 
tonian will be considered in the following. In this situa- 
tion, the unperturbed system has a conserved Hamilto- 
nian Ho, a Liouville operator iL^ = ^iB^d/dxi, with 
defined by Eq. 10), and an equilibrium distribution 
function given by pe — \\Z\\S{x)S{'Hq). The Liouville 
equation in the presence of the perturbation is 

^ = -iC^p-tLf{t)p, (33) 

where Q = iLq + k^- One can consider p = p^ + Sp 
and to linear order 

6p{t) = - f dTT{T)e-'^oit-r),L^p^ . (34) 
Jo 

The nonequilibrium average of Sb{t) — b{t) — {b)eq for any 
phase space variable is then 

6b{t) = / dTj=-{T)(l){t - t) , (35) 
Jo 

where 

(j)[t) = - j dxb{t)iL^pe (36) 

is the response function and b{t) — expliL^' t]b{0) since 
the compressibility Kc disappears when integrating by 
parts exp[i£^t] 0, uM- Now, in evaluating the action 
of iLp on Pe one can take full advantage of the fact that 
iLpS{x) — 0- Had one used iLp instead, spurious terms 
would have appeared. Thus 

iL^Pe = Pe[iL^ln\\Z\\-P{iL^no)] , (37) 



where it has been used the fact that iLp5{T-lo) « 
—PS (Ho) {iLp Ho) in the thermodynamic hmit [33l| . 
Hence, the response function for constrained systems 
takes the form 

cp{t)^f3{{iL^Ho-l3-'tL^\n\\Z\\)),, . (38) 

The correction factor arising from iL^ In ||.Z|| disappears 
if 

iLp = ^D,5/api (39) 

i 

but for general equations of motion [H IHil HI IH 
I23, 11^123 it must be considered. 



V. CONCLUSIONS 

The extension of Dirac's formalism allows one to treat 
correctly systems with holonomic constraints undergoing 
non-Hamiltonian dynamics. Non-Hamiltonian dynamics 
can be derived from generalized brackets or it can be 
more general and not be derivable from any brackets: 
both cases have been treated. Using generalized brack- 
ets, a Nose-Dirac phase space flow has been introduced 
and applied to derive the unbiasing factor when con- 
straints are used to sample rare events. It has been shown 
how to generalize Dirac's recipe when the dynamics can- 
not be obtained from brackets. Linear response theory 
of system with holonomic constraints subjected to gen- 
eral non-Hamiltonian perturbation has been illustrated. 
The use of Dirac's formalism makes spurious terms disap- 
pear from the response function. However, a correction 
coming from the measure of constrained phase space is 
present in general cases. Further work is required in order 
to assess the importance of this correction in numerical 
calculations on condensed matter systems. 

Equilibrium statistical mechanics and linear response 
theory of systems with nonholonomic constraints remain 
to be addressed. However, as a consequence of the anal- 
ysis presented in this paper, it can be suggested that a 
formalism, suitable for the linear response of such sys- 
tems, must project the spurious time evolution of the 
constrained variables out of both unperturbed and per- 
turbed dynamics. 
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APPENDIX A: NOSE-DIRAC FLOW AND RARE 
EVENTS 

Often one is interested in the calculation of conditional 
averages of some phase space function a(r) 



(a(r)} 



cond 



{a{T)Si^ir)~e))NVT 



The relation between conditional average IjAip and 
constrained averages HA2|I has been originally given in 
Ref. 0. In the present context, it is simply remarked 
that, since the formal manipulations are performed in 
the canonical ensemble, in order to be rigorous one needs 
the Nose-Dirac flow to have the correct distribution func- 
tion given in Eq. (EHJ- Having said that, one just needs 
(Al) the results of Ref 34], which show that / d^p| |Z| |(5(o-) cx 
1^11^/^, in order to re-write the constrained average as 



where (. . ■)nvt stands for an equilibrium average in the 
canonical ensemble. Molecular dynamics can be used to 
perform calculation where the condition cr(r) = ^ — = 
is treated as a holonomic constraint. However this au- 
tomatically brings in a constraints on the time variation 
^■(rjr) so that using the Nose-Dirac flow introduced be- 
fore one would get a constrained average, defined by 



(a(r))e, = 
where S{x) — S{cr)5{cr) 



{a{r)\\Z\\Six))NVT 

(II^II'^(X))WT 



(A2) 



(«(r))ct 



(a(r)||Z||i/2^(^))^^^ 



(A3) 



{\\Z\\^/^6icT))NVT ■ 
From this, one immediately obtains the result of Ref. 



(a(r)) 



cond 



(A4) 
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